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This paper presents an extension of a numerical algorithm for network flow analysis code 
to perform multi-dimensional flow calculation. The one dimensional momentum equation 
in network flow analysis code has been extended to include momentum transport due to 
shear stress and transverse component of velocity. Both laminar and turbulent flows are 
considered. Turbulence is represented by Prandtl’s mixing length hypothesis. Three 
classical examples (Poiseuille flow, Couette flow and shear driven flow in a rectangular 
cavity) are presented as benchmark for the verification of the numerical scheme. 


1.0 Introduction 

Traditionally, fluid network codes have been used for system level analyses, whereas 
Navier-Stokes(NS) codes have been used for component level analyses. Until recently, 
most attempts to merge the two methodologies have come from the NS side (i.e. using a 
NS code to perform a system level analysis). This approach brings the enormous 
overhead associated with such a code. The current approach, on the other hand, begins 
from the fluid network side. The system level code utilized in the current study is the 
Generalized Fluid System Simulation Program (GFSSP) [1]. 

The Generalized Fluid System Simulation Program was developed for the purpose of 
calculating pressure and flow distribution in a complex flow network associated with 
secondary flow in a liquid rocket engine turbopump. The code was developed to be a 
general purpose flow network solver so that generic networks could be modeled. A given 
fluid system is discretized into the nodes and branches. This practice is conceptually 
similar to the “staggered grid” practice of SIMPLE algoritm of Patankar & Spalding[2]. 
GFSSP employs a finite volume formulation of mass, momentum, and energy 
conservation equations in conjunction with the thermodynamic equations of state for real 
fluids. Mass, energy and specie conservation equations are solved at the nodes; the 
momentum conservation equations are solved in the branches. The system of equations 
describing the fluid network is solved by a hybrid numerical method that is a combination 
of the Newton-Raphson and successive substitution methods. Eighteen different 
resistance/source options are provided for modeling momentum sources or sinks in the 
branches. Two thermodynamic property programs, GASP-WASP[3,4] and GASPAK[5] 
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are integrated with the code to provide thermodynamic and thermophysical properties of 
real fluid. GFSSP’s system level capability has been extensively verified by comparing 
with test data [6-10]. 


2.0 Unstructured Finite Volume Grid 

The unstructured finite volume grid network for GFSSP is shown in Figure 1, which 
shows connectivity of five nodes with four branches. In this figure node-i is connected 
with four neighboring nodes (j = 1 to 4) . In structured coordinate systems the number of 
neighboring nodes are restricted to 2, 4 and 6 for one, two and three dimensional systems 
respectively. On the other hand for an unstructured system, there is no such restriction on 
the number of neighboring nodes. The index k represents fluid species. 



Fluid (k=2) 


Figure 1: Schematic of Nodes and Branches for an Unstructured Finite Volume Grid 


3.0 Conservation Equations 

3.1 Mass Conservation Equation 

The mass conservation equation for the i th node can be represented by: 


n 

S 

j = l 


m.. = 0 
ij 


( 1 ) 


Equation 1 implies that the net mass flow from a given node must equate to zero. In other 
words, the total mass flow rate into a node is equal to the total mass flow rate out of the 
node. 
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3.2 Momentum Conservation Equation 

The one dimensional form of momentum equation for every branch takes the following 
form: 
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Equation 2 represents the balance of fluid forces acting on a given branch. Inertia, 
pressure, gravity, friction and centrifugal forces are considered in the conservation 
equation. In addition to five forces, a source term S has been provided in the equation to 
input pump characteristics or to input power to pump in a given branch. If a pump is 
located in a given branch, all other forces except pressure are set to zero. The source 
term S is set to zero in all other cases. Figure 2 shows the schematic of a branch control 
volume. 



Figure 2: Schematic of a Branch Control Volume Showing the Gravity and Rotation 

Multi-dimensional conservation equations must account for the transport of mass, 
momentum and energy into and out of the control volume from all directions in space. 
Mass conservation equations (Equation 1) can account for such transport because each 
internal node can be connected with multiple neighboring nodes located in space in any 
arbitrary location (Figure 1). On the other hand, the momentum conservation equation 
(Equation 2) is one dimensional. Multi-dimensional momentum transport can be 
accounted for by incorporating two additional terms in the momentum equation. These 
terms include: a) momentum transport due to shear, and b) momentum transport due to 
the transverse component of velocity. 
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These two terms can be identified in the two-dimensional steady state Navier-Stokes 
equation which can be expressed as: 
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The first term on the left hand side of Equation 3 corresponds to the current inertia term 
described in Equation 2. The second term on the left hand side of Equation 3 
corresponds to the transverse momentum exchange. The first term on the right hand side 
corresponds to the pressure term in Equation 2. The second term on the right hand side 
of Equation 3 corresponds to the gravity term. The third term on the right hand side of 
the equation is negligible (based on an order of magnitude argument). The fourth term on 
the right hand side represents momentum transport due to shear. The next two sub- 
sections describe the implementation of shear and transverse momentum transport into 
the momentum equation of GFSSP. A more detailed description of the implementation 
of these terms for laminar flow is detailed by Schallhom [11]. 

3.3 Momentum Transport Due to Shear 


3.3.1 Laminar Flow 


Begin by examining the shear term (fourth term) of the Navier-Stokes Equation in more 
detail. First, consider the shear as a force instead of a force per unit volume by 
multiplying the volume by the shear term. 

p|^V = p-^ r (AxXAyXAz) = p-^-AxAz = p-^A Aear (4) 

dy 2 (Ay) Ay Ay 


Figure 3 represents a set of nodes and branches for which shear forces are exchanged. 
Let branch 12 represent the branch for which the shear force is to be calculated. 
Branches N12 and S12 represent the parallel branches which will be used to calculate the 
shear force on branch 12. Let YS be the distance between branches 12 and S12, and let 
YN be the distance between branches 12 and N12. Let AS be the shearing area between 
branches 12 and S 12, while AN is the shearing area between branches 12 and N12. 
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A differencing scheme that can account for non-orthogonality in node structure was used. 
Equation 5 represents the shear term for branch 12. The angle 0 represents the angle that 
adjacent branches make with respect to the referenced branch. 
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Equation 5 can be generalized to n-number of parallel branches at any position around 
branch 12 as shown in Equation 6. 
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( 6 ) 


Where the summation subscript i represents the i th parallel branch to branch 12. 

Next, the shear interaction with a neighboring parallel wall should be addressed. 
Suppose that adjacent to branch 12 is a wall that is approximately parallel to the branch. 
The angle O wa ii represents the angle between branch 12 and the wall. The wall has a 
velocity v^nd- The distance between the centerline of branch 12 and the wall is y wa ii and 
the shear area is A wa ». The expression for the shear effect of the wall on branch 12 is 
given in Equation 7. 


3u . 
p— A 
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If there were multiple walls adjacent to branch 12, then Equation 7 could be generalized 
into Equation 8. 



Branch 12-wall,^ 
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( 8 ) 
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Finally, combining the adjacent branch and the wall shear equations into one generalized 
equation for the i l " branch for which shear is to be calculated. Equation 9 represents the 
actual laminar shear formulation incorporated into GFSSP, where i is the current branch, 
npi is the number of parallel branches to branch i, and ns* is the number of parallel solid 
walls to branch i. 
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3.3.2 Turbulent Flow 


To model shear stress in a turbulent flow, a turbulence model must be employed. 
Separate methods must be used to model the interaction between adjacent branches and 
between branches and adjacent walls. The following two sections provide the details for 
the modeling of these interactions. 

3.3.2. 1 Branch-Branch Interaction 

Turbulent shear stress interaction between a branch and one of its neighboring branches 
is modeled using a modified form of the Mixing Length algebraic model proposed by 
Prandtl (see reference 5 for a description of Prandtl Mixing Length model). Referring to 
Figure 3 again, shear interaction between branch 12 and branch N12 is determined by 
their relative mean velocities, the shear area (AN) and the distance between the two 
branches (YN). In the original mixing length turbulence model proposed by Prandtl, 
viscosity (beginning with Equation 3) is replaced by an effective viscosity. This effective 
viscosity is defined as the sum of the viscosity and a “turbulent” viscosity: 


Ineffective — P + P 


turbulent * 


( 10 ) 


The turbulent viscosity, proposed by Prandtl, is defined by Equation 11, below. 


Peurb =pM 2 


du 

dy 


where, p = local density of the fluid, 

K = Prandtl’s mixing length constant (0.4), 
y = distance from the wall, 
u = local fluid velocity parallel to the wall. 


(ID 


Prandtl’s formulation requires knowledge of local positioning with respect to a wall(s); 
however, GFSSP’s formulation (outlined earlier) is fully unstructured. In order to 
implement the above approach into GFSSP, either additional information is required in 
the input file, or a modification to the definition of “y” in Equation 1 1 is needed. Based 
upon a desire to easily allow for an individual model to provide results for both laminar 
and turbulent approaches with minimal input change by the user, the latter approach of 
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modifying the definition of “y” was chosen. The new definition of “y” for Equation 1 1 is 
the distance between the branches, which is a required input for the laminar approach 
already. Therefore, for turbulent flow, Equation 4, 5, and 6 become Equations 12, 13, 
and 14. 
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Equation 14 represents the actual branch-branch turbulent shear formulation put into 
GFSSP. 


3. 3.2.2 Branch -Wall Interaction 

Turbulent shear stress interaction between a solid (wall) and an adjacent branch is 
modeled using the log law of the wall [12]. The log law of the wall utilizes a 
characteristic turbulence distance from the wall (y + ) and a characteristic velocity (u + ). y + 
is a function of wall shear stress, and u + is a function of y + and wall shear stress; therefore 
an iterative scheme is required to calculate wall shear stress. Equations 15, 16, and 17 
are the three equations that are iterated upon until a converged value for wall shear stress 
is achieved. 
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where, p = local fluid density 
p = local fluid viscosity 
y = local distance from wall 


(15) 


7 



(16) 


u 


+ 


y + 

- 3.05 + 5.o(log(y + )) 
5.5 + 2.5(log(y + )) 


if y + <5.0 
if 5.0 <y + <30.0 
if y + >30.0 


_n| 


JL) 

_P I 

V 



3u 

where, u = local velocity parallel to the wall 


(17) 


In order to initiate the iterative process, a seed value of the wall shear stress must be 
provided. The seed value of wall shear stress utilized is given by Equation 18: 
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initial guess 


pu 
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(18) 


Equations 15-18 represent the actual branch-wall turbulent shear stress formulation 
incorporated into GFSSP. 


3.4 Transverse Momentum Transport 

The transverse momentum component of Equation 3 can be expressed in terms of a force 
per unit volume. 

- 7 — V * ^ 77 — 7 ^- (A x)(A y)(A z) = (A x)(A z)(pvAu) = m^Au (19) 

dy (A y) 

Figure 4 represents a set of nodes and branches for which transverse momentum 
exchange will take place. Let the branch 12 represent the current branch which will 
receive transverse momentum from the surrounding branches. Branch S12 represents an 
adjacent parallel branch, while branches SI and S2 represent the adjacent normal 
branches. 



Figure 4: Branch and Node Schematic for Transverse Momentum Exchange 
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Now, examine the formulation for calculating the transverse momentum term for branch 
12. First, calculate the average mass flow rate for the adjacent normal branches: 


m s — + hi s2 ) 


( 20 ) 


Examining Figure 4, a positive transverse mass flow rate is defined as flow into the nodes 
corresponding to the branch in question. Based on this definition of transverse mass flow 
rate, calculate the transverse momentum term: 


dpvu v 

dy 


12 


(m S U,2 ^S^S12 ) *^s(**12 ^SI2 ) 


( 21 ) 


The parallel branch (S12) will contribute to the transverse momentum term of branch 12 
when m l2 >0 (since a positive transverse flow rate begins at S12 and ends at 12), and 
should have a negligible contribution when m, 2 <0 (i.e. transverse flow rate begins at 12 
and ends at SI 2). Equation 22, is the upwinding representation of the transverse 
momentum term for branch 12. 


V * u 12 (max||- ih s ,0)|)- u sl2 (max||0,m s ||) 
d y i2 


( 22 ) 


Equation 22 can be generalized for m-parallel branches around branch i, each with an 
angle By with respect to branch i, and n,j corresponding transverse connecting branches, 
each transverse branch with an angle of 0,jk with respect to branch i. Equation 23 
represents this generalized version of Equation 22. 
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Equation 23 represents the actual transverse momentum formulation put into GFSSP, 
where i is the current branch for which transverse momentum is being calculated, mi is 
the number of parallel branches which will be used to calculate transverse momentum, 
and n;j is the number of connecting transverse branches between the current branch i, and 
the j 01 parallel branch. 


4.0 Verification Results 


In order to verify proper implementation of the shear and transverse momentum 
components into GFSSP, three models were identified and developed. Two verification 
models were benchmarked for both laminar and turbulent flow. These two models are: 
two dimensional Poiseuille flow and two dimensional Couette flow. The third model was 
benchmarked with the laminar flow solution only. This third model is the two 
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dimensional shear driven flow in a square cavity. The following sections describe the 
models and presents the results. 

4.1 Poiseuille Flow Model 

Consider the flow between two fixed flat plates shown in Figure 5. The flow is pressure 
driven and assumed to be fully developed. The analytical solution for this situation can 
easily be derived for laminar flow. Figure 6 shows an approximate velocity profile for 
the laminar situation. 

A simple 3 node, 10 branch model (two sets of 5 parallel branches) was constructed to 
model the physical situation described above. The model is shown schematically in 
Figure 7. 
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top surface 


= 0 


upstream 


= 20 psi 
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bottom surface 


= 0 


^downstream P'* * 


Figure 5: Poiseuille Flow Physical Situation 



Figure 6: Poiseuille Flow Velocity Profile 


Top Wall Top Wall 



Bottom Wall Bottom Wall 


Figure 7: GFSSP Poiseuille Flow Model 
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The GFSSP’s flow through a restriction resistance option was used in the initial flow 
field calculation (for the first Newton-Raphson iteration, after which the shear will 
replace the friction factor calculation) for each of the branches. The bottom and top walls 
are fixed. 

Figure 8 shows a comparison between the velocity profiles for the known solutions 
(analytical for laminar and experimental for turbulent) and the GFSSP 3 node, 10 branch 
(5 parallel branch) model. As can be seen in Figure 8, the results of this crude GFSSP 
model compare vary favorably with the analytical solution. 



Figure 8: Poiseuille Flow Velocity Distribution 


4.2 Couette Flow 

Consider the flow between two fixed flat plates shown in Figure 9. The flow is shear 
driven and assumed to be fully developed. The solution to the laminar case is a linear 
velocity profile, illustrated in Figure 10. 
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Figure 9: Couette Flow Physical Situation 



Figure 10: Laminar Couette Flow Velocity Profile 

A simple 3 node, 10 branch model (two sets of 5 parallel branches) was constructed to 
model the physical situation. The model is shown schematically in Figure 11. 


Top Wall Top Wall 



Bottom Wall Bottom Wall 


Figure 1 1 : GFSSP Couette Flow Model 

As with the Poiseuille flow case, resistance option -02 was used in the initial flow field 
calculation (for the first Newton-Raphson iteration, after which the shear will replace the 
friction factor calculation) for each of the branches. The bottom walls are fixed, and the 
top walls are moving at a known velocity. 

Figure 12 shows the comparison between the velocity profiles for the known solutions 
(analytical for the laminar case, experimental results for the turbulent case) and the 
GFSSP 3 node, 10 branch (5 parallel branch) model. As can be seen in Figure 12, the 
results of this crude GFSSP model compare nearly identically with the analytical 
solution. 
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4.3 Shear Driven Flow in a Square Cavity 

Consider a square cavity as shown in Figure 13. The flow is induced by shear interaction 
at the top wall. The bottom and side walls are fixed. The top wall is moving to the right 
at constant speed. The corresponding Reynolds number for this situation is Re = 100. 

4.3.1 Benchmark Numerical Solution 

Due to the non-linearity of the governing differential equations, an analytical solution of 
this situation is not available. Instead of an analytical solution, a well known numerical 
solution by Odus Burggraf [13] was used as the benchmark. Burggraf used a 51x51 grid 
in his model of the square cavity. 

4.3.2 GFSSP Driven Cavity Model 

The GFSSP model of the driven cavity consists of 49 nodes (48 of which are internal) 
and 84 branches. 
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For numerical stability, one node (Node 1) was assigned to be a boundary node with an 
arbitrary pressure. A unit depth was assumed for the required areas. The model is shown 
schematically in Figure 14. As in the previous cases, GFSSP’s flow through a restriction 
resistance Option was used in the initial flow field calculation for all of the branches. 
The bottom and side walls are fixed. The top walls are moving to the right at known 
velocity. All parallel angles are 0°, and all transverse angles are 90°. 
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4.3.3 Results 


Figure 15 shows a comparison between the benchmark numerical solution and GFSSP 
7x7 node model velocity profiles along a vertical plane at the horizontal midpoint. As 
can be seen in Figure 15, the results of this crude GFSSP model compare very favorably 
with the benchmark numerical solution of Burggraf [13]. 



Figure 15: Shear Driven Square Cavity Centerline Velocity Distribution 


5.0 Conclusions 

This paper presents a numerical algorithm to extend a system level flow network code 
that was designed to solve one dimensional momentum equation to perform multi- 
dimensional flow calculation. The algorithm uses an identical mathematical framework 
for both system and component level analysis. The multi-dimensional features were 
incorporated by including additional momentum sources due to shear stress and transport 
of momentum due to transverse component of velocity. The simplicity and ease of the 
formulation can largely be attributed to the use of unstructured co-ordinate system. 
Excellent agreement with analytical solution was obtained for laminar flow in three 
benchmark problems: Poiseuille Flow, Couette flow and shear driven flow in a square 
cavity. Turbulence is modeled by an effective viscosity which is calculated from 
Prandtle’s mixing length theory. Numerical predictions compared well with turbulent 
Couette flow data. 
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